function [F, g, G] = DBV(x)
% DBV Discrete Boundry Value - 输出Discrete Boundry Value函数在x点的函数值，与梯度值。
% Call：
% [F, g, G] = DBV(x)
%    [F, g] = DBV(x)
%       [F] = DBV(x)
% 输入参数：
% x     ：给定的点，函数的n由x长度确定。
% 输出参数：
% F     ：函数值
% g     ：梯度值，不要求不计算
% G     ：Hessen矩阵，不要求不计算
    [stop, n] = check(x);
    if stop, error('x must be a real valued vector'), end
    if n <2, error('x must have length more than 1'),end
    x = x(:);
    
    r = zeros(n,1);
    h = 1/(n+1);
    r(1) = 2*x(1)-x(2) + h^2*(x(1)+h+1)^3/2;
    r(n) = 2*x(n)-x(n-1) + h^2*(x(1)+h*n+1)^3/2;
    for i = 2:n-1,
        t = h*i;
        r(i) = 2*x(i)-x(i-1) - x(i+1) + h^2*(x(1)+t+1)^3/2;
    end    
    F = sum(r.^2);
    % 计算梯度：    
    if nargout>1,
        pr = zeros(n,n);
        g  = zeros(n,1);
        pr(2,1) = -1;
        pr(1,1) = 2 + 3*h^2*(x(1)+h+1)/2; 
        pr(n-1,n) = -1;
        pr(n,n) = 2 + 3*h^2*(x(n)+h*n+1)/2; 
        for i=2:n-1,
            t  = h*i;
            pr(i-1,i) = -1;
            pr(i+1,i) = -1;
            pr(i,i) = 2 + 3*h^2*(x(i)+t+1)/2; 
        end
        for i =1:n,
            g = g + r(i)*pr(:,i);
        end
        g = 2*g;
    end
    % 计算Hessen矩阵
    if nargout>2,
        G = zeros(n, n);
        Hr = zeros(n,n,n);
        for i = 1:n,
            t = h*i;
            G(i,i) = G(i,i)+3*h^2*(x(i)+t+1)*r(i);
            for l = 1:n,
                for k = 1:n,
                    G(l,k) = G(l,k) + pr(l,i) * pr(k,i);
                end
            end
        end
        G = 2*G;
    end
end
%============  Auxiliary functions  ========================

function [err,n] = check(x)
% CHECK - check x
%   
    err = 0; sx = size(x); n = max(sx);
    if  (min(sx) ~= 1) | ~isreal(x) | any(isnan(x(:))) | isinf(norm(x(:))) 
        err = -1; 
    end    
end